Instala las librerías necesarias (copia y pega en la terminal; no descomentes la línea)…
# install.packages(
# c("tidyverse", "afex", "emmeans", "GGally", "effectsize", "performance", "see", "qqplotr")
# )
… y carga las librerías más usadas:
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✔ ggplot2 3.3.5 ✔ purrr 0.3.4
## ✔ tibble 3.1.7 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.0.1 ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
El modelado estadístico es difícil, por lo que en este notebook solo cubriremos algunos aspectos básicos. Si en el futuro te enfrentas a experimentos complejos, considera buscar ayuda.
To consult the statistician after an experiment is finished is often merely to ask him to conduct a post mortem examination. He can perhaps say what the experiment died of. (Ronald Fisher)
En el artículo clásico Statistical modeling: The two cultures encontramos
There are two cultures in the use of statistical modeling to reach conclusions from data. One assumes that the data are generated by a given stochastic data model. The other uses algorithmic models and treats the data mechanism as unknown. (Leo Breiman)
En este notebooks nos centraremos en modelos lineales y en los métodos de inferencia clásica. El siguiente notebook cubre modelos basados en aprendizaje automático.
(Otras lecturas interesante sobre las dos culturas es To explain or to predict?.)
El modelo básico sobre el que se construye gran parte de la “Estadística clásica” es el modelo de regresión lineal: \[y = a + b\cdot x + \epsilon\] donde \(\epsilon \sim \mathcal{N}(0, \sigma^2)\).
Es instructivo simular datos que sigan este modelo para entender el significado de la ecuación:
x = seq(-2, 2, 0.1)
expected_behaviour = 2 + 3 * x # a = 2 y b = 3
epsilon = rnorm(length(x), sd = 1)
y = expected_behaviour + epsilon
df = data.frame('data_x' = x, 'data_y' = y,
'expected' = expected_behaviour)
ggplot(df, aes(x = data_x)) +
geom_point(aes(x = data_x, y = data_y)) +
geom_line(aes(y = expected), col = 2)
# OR...
# plot(x, y)
# lines(x, expected_behaviour, col = 2)
Podemos pensar que expected_behaviour (\(a + b\cdot x\)) es el comportamiento medio
o esperado en la población de interés, mientras que epsilon
(\(\epsilon\)) representan
fluctuaciones aleatorias (bien debidas a errores del proceso de medición
o que el proceso estudiado tiene cierta aleatoriedad). Además, la
relación entre \(x\) e \(y\) es muy específica. Si aumenta
(disminuye) \(x\) aumenta (disminuye)
\(y\). Además, un incremento de una
unidad en \(x\) siempre produce el
mismo incremento en \(y\) (ídem si
\(x\) disminuye). Se dice que la
relación entre \(x\) e \(y\) es lineal.
La primera pregunta a la que nos enfrentamos es la siguiente.
Dados los datos \((x, y)\),
¿podemos estimar expected behaviour?
# 1) crear un modelo lineal
naive_model = lm(data_y ~ data_x, data = df)
# 2) obtener estimaciones de a y b
summary(naive_model)
##
## Call:
## lm(formula = data_y ~ data_x, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.79042 -0.61648 -0.07474 0.78419 1.93933
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.0753 0.1439 14.42 <2e-16 ***
## data_x 2.8666 0.1216 23.57 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9215 on 39 degrees of freedom
## Multiple R-squared: 0.9344, Adjusted R-squared: 0.9327
## F-statistic: 555.5 on 1 and 39 DF, p-value: < 2.2e-16
# 3) Obtener predicciones del modelo lineal
preds = predict(naive_model, interval = "confidence")
# 4) visualizar el ajuste
data_and_preds = df %>% bind_cols(preds)
ggplot(data_and_preds, aes(x = data_x)) +
geom_point(aes(y = data_y)) +
geom_line(aes(y = expected, col = "Expected")) +
geom_line(aes(y = fit, col = "Predicted")) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.2, fill = "black")
Dado que las estimaciones incorporan el error de estimación, es posible realizar inferencia acerca de la significación de los parámetros.
Un estudiante de biología desea determinar la relación entre temperatura ambiente y frecuencia cardíaca en la rana leopardo, Rana pipiens. Para ello, manipula la temperatura en incrementos de 2ºC que van desde 2ºC a 18ºC, registrando la frecuencia cardíaca (pulsaciones por minuto) en cada intervalo. Los datos están disponibles en “hr.csv”.
# 1) leer los datos
frog_df = read.table("data/hr.csv", header = TRUE)
# 2) Crear un modelo lineal
frog_model = lm(heart_rate ~ temperature, data = frog_df)
# 3) Inferencia
summary(frog_model)
##
## Call:
## lm(formula = heart_rate ~ temperature, data = frog_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3389 -1.7889 -0.6889 1.7611 5.0111
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1389 1.9170 1.116 0.301
## temperature 1.7750 0.1703 10.421 1.63e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.639 on 7 degrees of freedom
## Multiple R-squared: 0.9394, Adjusted R-squared: 0.9308
## F-statistic: 108.6 on 1 and 7 DF, p-value: 1.629e-05
El diseño experimental y los resultados de la inferencia en el ejemplo de las ranas nos invitan a concluirque “el aumento de la temperatura causa un incremento de la frecuencia cardíaca”. Sin embargo, esto no es correcto. Por muy fuerte que parezca la relación entre las variables \(x\) e \(y\), nunca debemos interpretar una variable como la causa de la otra. Una relación significativa entre \(x\) e \(y\) puede ocurrir por varios motivos:
Por otra parte, respecto a los p-valores…
There is some debate among statisticians and researchers about the appropriateness of P values, and that the term “statistical significance” can be misleading. If you have a small P value, it only means that the effect being tested is unlikely to be explained by chance variation alone, in the context of the current study and the current statistical model underlying the test. If you have a large P value, it only means that the observed effect could plausibly be due to chance alone: it is wrong to conclude that there is no effect (emmeans package authors)
En general, hay consenso en que debe abandonarse la interpretación dicotómica del p-valor (efecto significativo Vs no-significativo, sobre todo teniendo en cuenta que se basan en un threshold arbitrario) y que debe favorecerse los resultados basados en intervalos de confianza y tamaños de los efectos (¡incluso si no son significativos!)
For example, a study on the effects of two different ambient temperatures on paramecium diameter returning an effect size of 20 µm and a p-value of 0.1, if centred on p-value interpretation would conclude ‘no effect’ of temperature, despite the best supported effect size being 20, not 0. An interpretation based on effect size and confidence intervals could, for example, state: ‘Our results suggest that paramecium kept at the lower temperature will be on average 20 µm larger in size, however a difference in size ranging between −4 and 50 µm is also reasonably likely’. (…), the latter approach acknowledges the uncertainty in the estimated effect size while also ensuring that you do not make a false claim either of no effect if p > 0.05, or an overly confident claim. (Lewis G. Halsey, The reign of p-value is over)
summary(frog_model)
##
## Call:
## lm(formula = heart_rate ~ temperature, data = frog_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3389 -1.7889 -0.6889 1.7611 5.0111
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1389 1.9170 1.116 0.301
## temperature 1.7750 0.1703 10.421 1.63e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.639 on 7 degrees of freedom
## Multiple R-squared: 0.9394, Adjusted R-squared: 0.9308
## F-statistic: 108.6 on 1 and 7 DF, p-value: 1.629e-05
confint(frog_model)
## 2.5 % 97.5 %
## (Intercept) -2.394015 6.671792
## temperature 1.372241 2.177759
Los resultados sugieren que el ritmo_cardiaco de las ranas se incrementará 1.77 bpms por cada grado Celsius, si bien un aumento de la frecuencia cardiaca en el rango (1.37, 2.17) es igualmente verosímil.
Evidentemente cualquier interpretación está supeditada a que el modelo sea correcto. Debemos ser muy cuidadosos a la hora de verificar que se cumplan las asunciones del modelo de regresión lineal. Podemos usar el acrónimo LINE para recordar las asunciones más importantes del modelo: Linear, Independent, Normal, Equal variances.
naive_modelplot(naive_model, ask = FALSE)
# Comprobar la normalidad con qqplot puede ser difícil. Podemos apoyarnos en
# performance::check_normality
library("performance")
# check_normality corre shapiro.test, pero tal y como resalta la documentación
# "this formal test almost always yields significant results for the distribution
# of residuals and visual inspection (e.g. Q-Q plots) are preferable."
is_norm = check_normality(naive_model)
# Para hacer la inspección visual
plot(is_norm)
plot(is_norm, type = "qq")
plot(is_norm, type = "qq", detrend=TRUE)
Després et al. 1 señalan que el tejido adiposo (AT) se asocia con complicaciones metabólicas consideradas como factores de riesgo para la enfermedad cardiovascular. Es importante, afirman, medir la cantidad de AT intraabdominal como parte de la evaluación del riesgo de enfermedad cardiovascular de un individuo. Després y sus colegas realizaron un estudio para predecir la cantidad de AT abdominal a partir de simples mediciones antropométricas.
Entre las medidas tomadas en cada sujeto, se incluyó la AT abdominal obtenida por TC y la circunferencia de la cintura (datos en Wat.csv’’). Una pregunta de interés es cómo de bien se puede predecir y estimar la AT abdominal a partir del conocimiento de la circunferencia de la cintura. Construye un modelo lineal y plotea sus predicciones. ¿Existe alguna relación entre AT y la circunferencia de la cintura?
at_df = read.table("data/at.csv", header=TRUE, sep = ",")
names(at_df) = c("subject", "waist", "AT")
at_model = lm(AT ~ waist, data = at_df)
at_with_preds =
at_df %>%
mutate(fit = predict(at_model))
confint(at_model)
## 2.5 % 97.5 %
## (Intercept) -259.190053 -172.77292
## waist 2.993689 3.92403
ggplot(at_with_preds, aes(x = waist, y = AT)) +
geom_point() +
geom_line(aes(y = fit, col = "Predictions")) +
xlab("waist circumference (cm)") +
ylab("deep abdominal AT (cm2)")
¿Se cumplan las asunciones de la regresión lineal en el problema de AT abdominal?
plot(at_model, ask = FALSE)
# No, el ruido es heterocedastico (no es el mismo para todos los valores de x):
# ver plot 3.
Para lidiar con situaciones como la ilustrada en el gráfico de “correlation is not causation” (tiburones Vs helados) necesitamos emplear modelos de regresión múltiple, dado que estos permiten “controlar” las variables de confusión. Crear un modelo de regresión múltiple es análogo al caso unidimensional…
Jansen y Keller[^2] utilizaron la edad y el nivel de educación para predecir la capacidad de dirigir la atención (CDA) en sujetos ancianos. CDA se refiere a los mecanismos neuronales que enfocan la mente en lo que es significativo al tiempo que bloquea las distracciones. El estudio recopiló información sobre 71 mujeres mayores con estado mental normal. En este estudio las puntuaciones más altas se corresponden con un mejor funcionamiento de la atención. Las mediciones en CDA, edad en años y el nivel de educación (años de escolaridad) están recogidos en “cda.csv”. Crea un modelo de regresión múltiple para predecir la puntuación de CDA.
[^2] D. Jansen, M. Keller, Çognitive Function in Community-Dwelling Elderly Women”, Journal of Gerontological Nursing, 29 (2003), 34-43.
library("GGally")
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
df = read.csv("data/cda.csv")
ggpairs(df)
cda_model = lm(CDA ~ AGE + EDLEVEL, df)
# cda_model = lm(CDA ~ ., df) # equivalente al comando anterior
df$predictions = predict(cda_model)
summary(cda_model)
##
## Call:
## lm(formula = CDA ~ AGE + EDLEVEL, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.9804 -2.2125 -0.0761 2.2824 9.1230
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.49407 4.44297 1.237 0.220498
## AGE -0.18412 0.04851 -3.795 0.000316 ***
## EDLEVEL 0.61078 0.13565 4.503 2.7e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.134 on 68 degrees of freedom
## Multiple R-squared: 0.3706, Adjusted R-squared: 0.3521
## F-statistic: 20.02 on 2 and 68 DF, p-value: 1.454e-07
plot(cda_model, ask=FALSE)
# todo parece correcto!
Hasta ahora solo hemos usado datos continuos, pero nada evita usar datos categóricos como predictores. ¡Ojo! Los coeficientes asociados a datos categóricos no deben interpretarse como una pendiente.
Construye un modelo de regresión lineal para predecir el peso de una persona a partir de los datos contenidos en “antrop.csv”. Interpreta los coeficientes de la regresión.
antrop = read.csv("data/antrop.csv")
antrop = mutate(antrop, male = factor(male))
antrop_model = lm(weight ~ height + male, data = antrop)
summary(antrop_model)
##
## Call:
## lm(formula = weight ~ height + male, data = antrop)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.5982 -3.5488 0.3738 3.2563 16.3890
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -29.46361 7.31841 -4.026 6.89e-05 ***
## height 0.46974 0.04886 9.615 < 2e-16 ***
## male1 1.23216 0.74242 1.660 0.0978 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.353 on 368 degrees of freedom
## Multiple R-squared: 0.36, Adjusted R-squared: 0.3565
## F-statistic: 103.5 on 2 and 368 DF, p-value: < 2.2e-16
antrop_preds = antrop %>% mutate(fit = predict(antrop_model))
ggplot(antrop_preds, aes(x = height, col=male)) +
geom_point(aes(y = weight)) +
geom_line(aes(y = fit), lwd = 3)
summary(antrop_model)
##
## Call:
## lm(formula = weight ~ height + male, data = antrop)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.5982 -3.5488 0.3738 3.2563 16.3890
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -29.46361 7.31841 -4.026 6.89e-05 ***
## height 0.46974 0.04886 9.615 < 2e-16 ***
## male1 1.23216 0.74242 1.660 0.0978 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.353 on 368 degrees of freedom
## Multiple R-squared: 0.36, Adjusted R-squared: 0.3565
## F-statistic: 103.5 on 2 and 368 DF, p-value: < 2.2e-16
El modelo se puede escribir como \[weight = -29 + 0.47 * height + 1.23 * is-male\] por lo que 1.23 significa que, de media y para una misma altura, los hombres pesan 1.23 Kg más que las mujeres (hemos ajustado por el efecto de la altura).
Usa intervalos de confianza para interpretar los resultados de la regresión.
confint(antrop_model)
## 2.5 % 97.5 %
## (Intercept) -43.8547524 -15.0724663
## height 0.3736688 0.5658194
## male1 -0.2277644 2.6920896
Aunque es cierto que de media los hombres pesan 1.23 Kg que las mujeres, ¡no podemos descartar que las mujeres pesen más que los hombres! (hasta 0.22 Kg).
Los datos contenidos en “howell1.csv” son datos censales parciales del área !Kung San compilados a partir de entrevistas realizadas a finales de la década de 1960. Crea un modelo para predecir el peso de los individuos a partir de la altura y el sexo. Evalúa la bondad del modelo.
howell = read.csv("data/howell1.csv", sep = ";")
howell = mutate(howell, male = factor(male))
howell_model = lm(weight ~ height + male, data = howell)
summary(howell_model)
##
## Call:
## lm(formula = weight ~ height + male, data = howell)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.4420 -3.6612 0.0892 3.5143 14.2151
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -33.847168 1.093840 -30.943 <2e-16 ***
## height 0.499848 0.007825 63.875 <2e-16 ***
## male1 0.734532 0.432259 1.699 0.0898 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.984 on 541 degrees of freedom
## Multiple R-squared: 0.8858, Adjusted R-squared: 0.8853
## F-statistic: 2097 on 2 and 541 DF, p-value: < 2.2e-16
howell_preds = howell %>% mutate(fit = predict(howell_model))
ggplot(howell_preds, aes(x = height, col=male)) +
geom_point(aes(y = weight)) +
geom_line(aes(y = fit), lwd = 3)
Sin ni siquiera usar plot(howell_model) ya somos capaces
de ver que el ajuste es malo… cualquier conclusión basada en un modelo
erróneo será errónea (garbage in, garbage out).
El dataset iris (puedes obtenerlo con
data(iris)) contiene medidas del sépalo y pétalo de varias
especies de iris. Construye un modelo lineal para predecir la longitud
del sépalo únicamente en función de la especie. Interpreta los
coeficientes de la regresión.
data(iris)
iris_model = lm(Sepal.Length ~ Species, iris)
print(summary(iris_model))
##
## Call:
## lm(formula = Sepal.Length ~ Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6880 -0.3285 -0.0060 0.3120 1.3120
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0060 0.0728 68.762 < 2e-16 ***
## Speciesversicolor 0.9300 0.1030 9.033 8.77e-16 ***
## Speciesvirginica 1.5820 0.1030 15.366 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5148 on 147 degrees of freedom
## Multiple R-squared: 0.6187, Adjusted R-squared: 0.6135
## F-statistic: 119.3 on 2 and 147 DF, p-value: < 2.2e-16
iris_preds = iris %>% mutate(fit = predict(iris_model))
ggplot(iris_preds, aes(x=Species, fill = Species)) +
geom_boxplot(aes(y=Sepal.Length)) +
geom_point(aes(y = fit), shape=4, size=3)
Al interpretar los coeficientes de la regresión, observamos que
lm ha tomado como referencia la especie
setosa. Esto se puede observar usando
contrasts.
contrasts(iris$Species)
## versicolor virginica
## setosa 0 0
## versicolor 1 0
## virginica 0 1
Es decir el modelo es \[sepal = mean-setosa-sepal + 0.93 * is-versicolor + 1.58 * is-virginica.\]
Sin embargo, podríamos reescribir el modelo de otra forma de forma que los coeficientes tengan otro significado. Un ejemplo sencillo sería: \[sepal = mean-versicolor-sepal + \alpha_1 * is-setosa + \alpha_2 * is-virginica.\]
En este caso, simplimente estamos variando la especie de referencia.
De hecho, contrast se puede modificar para usar como
referencia otro nivel del factor:
contrasts(iris$Species) = cbind("_is_setosa" = c(1, 0, 0),
"_is_virginica" = c(0, 0, 1))
iris_model_2 = lm(Sepal.Length ~ Species, iris)
print(summary(iris_model_2))
##
## Call:
## lm(formula = Sepal.Length ~ Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6880 -0.3285 -0.0060 0.3120 1.3120
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.9360 0.0728 81.536 < 2e-16 ***
## Species_is_setosa -0.9300 0.1030 -9.033 8.77e-16 ***
## Species_is_virginica 0.6520 0.1030 6.333 2.77e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5148 on 147 degrees of freedom
## Multiple R-squared: 0.6187, Adjusted R-squared: 0.6135
## F-statistic: 119.3 on 2 and 147 DF, p-value: < 2.2e-16
Lo interesante es que podemos ajustar los contrastes de forma que respondan a nuestras preguntas científicas. En general, estos contrastes deben ser ortogonales.
Imagínemonos el siguiente universo paralelo. En este universo paralelo solo existe la especie setosa. Una empresa de ingeniería genética te contrata para crear nuevas especies con un sépalo más grande. Desarrollas un método conocido como “Método V”, que tiene dos variantes “V-I” y “V-II”. Los experimentos con estas variantes dan lugar a dos nuevas especies que llamas versicolor (V-I) y virginica (V-II). Te planteas dos preguntas científicas:
levels(iris$Species)
## [1] "setosa" "versicolor" "virginica"
contrasts(iris$Species)
## _is_setosa _is_virginica
## setosa 1 0
## versicolor 0 0
## virginica 0 1
# Usamos versicolor como referencia
contrasts(iris$Species) = cbind(
"V - setosa" = c(-2, 1, 1),
"I - II" = c(0, 1, -1)
)
contrasts(iris$Species)
## V - setosa I - II
## setosa -2 0
## versicolor 1 1
## virginica 1 -1
v_model = lm(Sepal.Length ~ Species, iris)
summary(v_model)
##
## Call:
## lm(formula = Sepal.Length ~ Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6880 -0.3285 -0.0060 0.3120 1.3120
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.84333 0.04203 139.020 < 2e-16 ***
## SpeciesV - setosa 0.41867 0.02972 14.086 < 2e-16 ***
## SpeciesI - II -0.32600 0.05148 -6.333 2.77e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5148 on 147 degrees of freedom
## Multiple R-squared: 0.6187, Adjusted R-squared: 0.6135
## F-statistic: 119.3 on 2 and 147 DF, p-value: < 2.2e-16
confint(v_model)
## 2.5 % 97.5 %
## (Intercept) 5.7602675 5.9263991
## SpeciesV - setosa 0.3599303 0.4774031
## SpeciesI - II -0.4277344 -0.2242656
El método V parace producir sépalos más grandes. Por otra parte, V-II es mejor que V-I.
LLegados a este punto, ¡ya hemos cubierto el 90% de los contenidos habituales de un curso de estadística habitual! Aunque parezca mentira, ya hemos hecho análisis tan complejos como
antrop.csv o howell,Desde una perspectiva moderna, todos los análisis clásicos (T-test, anova, ancova) pueden considerarse como simples modelos de regresión. Esto demuestra el poder unificador de esta perspectiva. En las siguientes secciones revisaremos sin embargo estos modelos clásicos para afianzar la conexión con los modelos de regresión y profundizar en algunas cosas que nos han quedado en el tintero.
¿Depende la altura de los !Kung adultos del sexo del inviduo? Emplea los datos en “howell1.csv”.
adult_howell =
howell %>% filter(age >= 18) %>%
mutate(male = factor(male)) %>%
mutate(sex = fct_recode(male, "male" = "1", "female" = "0"))
height_model = lm(height ~ sex, data = adult_howell)
summary(height_model)
##
## Call:
## lm(formula = height ~ sex, data = adult_howell)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.6585 -3.4635 0.2965 3.4715 18.7115
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 149.5135 0.4049 369.25 <2e-16 ***
## sexmale 10.8450 0.5914 18.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.537 on 350 degrees of freedom
## Multiple R-squared: 0.49, Adjusted R-squared: 0.4885
## F-statistic: 336.3 on 1 and 350 DF, p-value: < 2.2e-16
Apoya tus conclusiones dibujando el histograma de las alturas por sexo
ggplot(adult_howell, aes(x = height, fill=sex)) + geom_density(alpha=0.2)
¿Cumple el modelo las asunciones de la regresión lineal?
hist(resid(height_model))
Lo cierto es que, aunque correcto, nuestra aproximación no es la
habitual a la hora de comparar la media de dos
poblaciones. Si las poblaciones son normales (o si podemos
invocar el Teorema Central del Límite), podemos hacer uso del
t.test.
t.test(height ~ sex, adult_howell)
##
## Welch Two Sample t-test
##
## data: height by sex
## t = -18.148, df = 323, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
## -12.020596 -9.669319
## sample estimates:
## mean in group female mean in group male
## 149.5135 160.3585
El fichero “hc.csv” contiene los datos de un experimento para verificar si una hormana de crecimiento funciona o no. ¿Existen diferencias significativas entre el grupo de control y el de tratamiento?
library("effectsize")
## Registered S3 method overwritten by 'parameters':
## method from
## format.parameters_distribution datawizard
hc = read.csv("data/hc.csv")
t.test(height_inc_cm ~ group, hc, alternative="less")
##
## Welch Two Sample t-test
##
## data: height_inc_cm by group
## t = -7.6007, df = 19998, p-value = 1.537e-14
## alternative hypothesis: true difference in means between group control and group treatment is less than 0
## 95 percent confidence interval:
## -Inf -0.0836217
## sample estimates:
## mean in group control mean in group treatment
## 5.989647 6.096364
cohens_d(height_inc_cm ~ group, data = hc, alternative = "less")
## Cohen's d | 95% CI
## -------------------------
## -0.11 | [-Inf, -0.08]
##
## - Estimated using pooled SD.
## - One-sided CIs: lower bound fixed at (-Inf).
# o bien...
height_test = t.test(
hc %>% filter(group == "control") %>% pull(height_inc_cm),
hc %>% filter(group == "treatment") %>% pull(height_inc_cm),
alternative = "less"
)
cohens_d(height_test)
## Cohen's d | 95% CI
## -------------------------
## -0.11 | [-Inf, -0.08]
##
## - Estimated using un-pooled SD.
## - One-sided CIs: lower bound fixed at (-Inf).
# ¡o bien podemos dejar que la librería elija el tamaño del efecto adecuado!
effectsize(height_test)
## Cohen's d | 95% CI
## -------------------------
## -0.11 | [-Inf, -0.08]
##
## - Estimated using un-pooled SD.
## - One-sided CIs: lower bound fixed at (-Inf).
Siempre debemos considerar el tamaño del efecto, bien en base a nuestro conocimiento o bien usando los estadísticos adecuados como Cohen’s D. Para Cohen’s D la heurística es:
En los datos de adult_howell, ¿podemos concluir que los
hombres son más altos que las mujeres? ¿Cuál es el tamaño del
efecto?
t.test(height ~ sex, adult_howell, alternative="less")
##
## Welch Two Sample t-test
##
## data: height by sex
## t = -18.148, df = 323, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group female and group male is less than 0
## 95 percent confidence interval:
## -Inf -9.8592
## sample estimates:
## mean in group female mean in group male
## 149.5135 160.3585
cohens_d(height ~ sex, data = adult_howell, alternative = "less")
## Cohen's d | 95% CI
## -------------------------
## -1.96 | [-Inf, -1.74]
##
## - Estimated using pooled SD.
## - One-sided CIs: lower bound fixed at (-Inf).
John M. Morton et al. [^3] examinaron la función de la vesı́cula biliar antes y después de la fundoplicatura, una cirugı́a para detener el reflujo. Los autores midieron la funcionalidad de la vesı́cula biliar calculando la fracción de eyección de la vesı́cula biliar (GBEF) antes y después de la fundoplicatura. El objetivo de la fundoplicatura es aumentar la GBEF, que se mide como un porcentaje. ¿Hay evidencia para concluir que la fundoplicatura aumenta el funcionamiento de la GBEF? Datos en “gbef_long.txt” (o “gbef.txt”, para un reto).
# Los datos en gbef.txt no están en el formato adecuado para un análisis estadístico...
# ¡hay que luchar para ordenarlos de forma correcta!
gbef_df = read.table("data/gbef.txt")
gbef_df[3, ] = c("ID", "", 1:(ncol(gbef_df) - 2))
gbef_df = t(gbef_df)
gbef_df = as.data.frame(gbef_df[-(1:2), ]) %>% setNames(gbef_df[1, ])
gbef_df =
gbef_df %>% mutate(
Preop = as.numeric(Preop),
Postop = as.numeric(Postop),
ID = factor(ID)
)
gbef_df_long =
gbef_df %>%
pivot_longer(Preop:Postop, names_to="class", values_to = "gbef")
gbef_df_long = read.table("data/gbef_long.txt", header = TRUE)
gbef_df_long =
gbef_df_long %>%
arrange(ID)
head(gbef_df_long)
## ID class gbef
## 1 1 Preop 22.0
## 2 1 Postop 63.5
## 3 2 Preop 63.3
## 4 2 Postop 91.5
## 5 3 Preop 96.0
## 6 3 Postop 59.0
preop = gbef_df_long %>% filter(class == "Preop") %>% pull(gbef)
postop = gbef_df_long %>% filter(class == "Postop") %>% pull(gbef)
gbef_test = t.test(postop, preop, paired = TRUE, alternative = "greater")
print(gbef_test)
##
## Paired t-test
##
## data: postop and preop
## t = 1.9159, df = 11, p-value = 0.04086
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
## 1.131919 Inf
## sample estimates:
## mean of the differences
## 18.075
cohens_d(gbef_test)
## Cohen's d | 95% CI
## -----------------------
## 0.55 | [0.03, Inf]
##
## - One-sided CIs: upper bound fixed at (Inf).
diff = postop - preop
hist(diff)
ANOVA permite comparar más de dos grupos entre sí (¡como ya hicimos con iris!). Asumiendo los peligros de dar recetas generales, en una primera aproximación podemos seguir los siguientes pasos:
Fíjate que ya hemos cubierto casi todos los pasos en los ejemplos de regresión. Los pasos novedosos son:
Repitamos el análisis realizado para el método V, incorporando además los nuevos pasos.
library("car") # Anova
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
data("iris")
# 1) Visualizar
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
ggplot(iris, aes(x=Species, y = Sepal.Length, fill=Species)) + geom_boxplot() +
coord_flip()
# 2) Especificar contrastes si tenemos alguna hipótesis específica
# (no te preocupes de los denominadores)
contrasts(iris$Species) = cbind(
"V - setosa" = c(-2, 1, 1) / 3,
"I - II" = c(0, 1, -1) / 2
)
contrasts(iris$Species)
## V - setosa I - II
## setosa -0.6666667 0.0
## versicolor 0.3333333 0.5
## virginica 0.3333333 -0.5
v_model = lm(Sepal.Length ~ Species, iris)
# 3) Correr ANOVA: test omnibus
v_model_aov = Anova(v_model, type = 3)
En el caso de modelos ANOVA, existen varios tamaños de efecto. Entre los mas conocidos están eta-squared y omega-squared. Para complicar aún más las cosas, las heurísticas para clasificar el tamaño como pequeño, mediano o largo son distintas que para Cohen’s d. En el caso que nos ocupa, para eta-squared o omega-squared tendríamos:
En general, consulta el siguiente FAQ para consultar las heurísticas.
eta_squared(v_model_aov) # o effectsize(v_model_aov)
## For one-way between subjects designs, partial eta squared is equivalent to eta squared.
## Returning eta squared.
## # Effect Size for ANOVA
##
## Parameter | Eta2 | 95% CI
## -------------------------------
## Species | 0.62 | [0.54, 1.00]
##
## - One-sided CIs: upper bound fixed at (1).
omega_squared(v_model_aov) # omega-squared se supone que está menos sesgado que eta_squared
## For one-way between subjects designs, partial omega squared is equivalent to omega squared.
## Returning omega squared.
## # Effect Size for ANOVA
##
## Parameter | Omega2 | 95% CI
## ---------------------------------
## Species | 0.61 | [0.53, 1.00]
##
## - One-sided CIs: upper bound fixed at (1).
Antes de seguir adelante deberíamos comprobar que se cumplan las asunciones de ANOVA. Dejémoslo para otro ejemplo y sigamos adelante.
En el código anterior, ANOVA nos indica que alguna(s) de las especies tiene(n) un sépalo distinto al de lasotras especies… ¡Sin embargo no nos dice cuáles son estas especies!
Para descubrirlo podemos usar contrastes …
# 4.a)
summary(v_model)
##
## Call:
## lm(formula = Sepal.Length ~ Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6880 -0.3285 -0.0060 0.3120 1.3120
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.84333 0.04203 139.020 < 2e-16 ***
## SpeciesV - setosa 1.25600 0.08916 14.086 < 2e-16 ***
## SpeciesI - II -0.65200 0.10296 -6.333 2.77e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5148 on 147 degrees of freedom
## Multiple R-squared: 0.6187, Adjusted R-squared: 0.6135
## F-statistic: 119.3 on 2 and 147 DF, p-value: < 2.2e-16
confint(v_model)
## 2.5 % 97.5 %
## (Intercept) 5.7602675 5.9263991
## SpeciesV - setosa 1.0797908 1.4322092
## SpeciesI - II -0.8554688 -0.4485312
… o usar tests post-hoc. La idea básica de los tests post-hoc es fácil de entender: dado que sé que existe alguna diferencia entre las especies, voy a comparar todas las especies entre sí. El problema es que esto dispara el error tipo I muy rápidamente, por lo que tenemos que ser más conservadores a la hora de aceptar una diferencia como significativa. Esto lo logramos con distintos métodos de ajuste de p-valores. Fíjate que el test omnibus sirve como una primera barrera protectora antes de lanzarnos a hacer comparaciones múltiples.
Entre los métodos de ajuste, podemos distinguir entre 1. Métodos centrados en controlar el family-wise error rate (FWER), cuyo credo es “Si cometo un solo error tipo I todas mis conclusiones se desmoronarán”. 2. Métodos centrados en controlar el false discovery rate (FDR), que se corresponde con el credo (considerablemente más optimista) “vamos a intentar estimar cuántos errores tipo I cometo y a no pasarme de cierto número (pero no pasa nada si hay algún error)”.
Podemos emplear R base para realizar los tests
post-hoc…
# Bonferroni es bastante conservador, pero es un ajuste muy conocido
pairwise.t.test(iris$Sepal.Length, iris$Species, p.adjust.method = "bonferroni")
##
## Pairwise comparisons using t tests with pooled SD
##
## data: iris$Sepal.Length and iris$Species
##
## setosa versicolor
## versicolor 2.6e-15 -
## virginica < 2e-16 8.3e-09
##
## P value adjustment method: bonferroni
# Los métodos fdr son "BH" (aka "fdr") and "BY".
pairwise.t.test(iris$Sepal.Length, iris$Species, p.adjust.method = "fdr")
##
## Pairwise comparisons using t tests with pooled SD
##
## data: iris$Sepal.Length and iris$Species
##
## setosa versicolor
## versicolor 1.3e-15 -
## virginica < 2e-16 2.8e-09
##
## P value adjustment method: fdr
… o bien el paquete emmeans (que tiene ciertas ventajas,
como veremos):
library("emmeans")
##
## Attaching package: 'emmeans'
## The following object is masked from 'package:GGally':
##
## pigs
v_model_emms = emmeans(v_model, "Species")
pairs(v_model_emms, adjust="bonferroni", infer=c(TRUE, TRUE))
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## setosa - versicolor -0.930 0.103 147 -1.179 -0.681 -9.033 <.0001
## setosa - virginica -1.582 0.103 147 -1.831 -1.333 -15.366 <.0001
## versicolor - virginica -0.652 0.103 147 -0.901 -0.403 -6.333 <.0001
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: bonferroni method for 3 tests
pairs(v_model_emms, adjust="fdr", infer = c(TRUE, TRUE))
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## setosa - versicolor -0.930 0.103 147 -1.179 -0.681 -9.033 <.0001
## setosa - virginica -1.582 0.103 147 -1.831 -1.333 -15.366 <.0001
## versicolor - virginica -0.652 0.103 147 -0.901 -0.403 -6.333 <.0001
##
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 3 estimates
## P value adjustment: fdr method for 3 tests
# una de las ventajas de emmeans es que podemos calcular tamaños de efecto
# para las medias aunque, tal y como se menciona en la documentación:
# "there is substantial disagreement among practitioners on what is the appropriate
# sigma to use in computing effect sizes; or, indeed, whether any effect-size
# measure is appropriate for some situations"
eff_size(v_model_emms, sigma=sigma(v_model), edf = 147)
## contrast effect.size SE df lower.CL upper.CL
## setosa - versicolor -1.81 0.226 147 -2.25 -1.360
## setosa - virginica -3.07 0.269 147 -3.60 -2.542
## versicolor - virginica -1.27 0.213 147 -1.69 -0.845
##
## sigma used for effect sizes: 0.5148
## Confidence level used: 0.95
Una desventaja de codificar los contrastes con contrasts
es que hay que prestar atención a los detalles (y el diablo esta en los
detalles). Por ejemplo, ¿de dónde sale el 3 del contraste
"V - setosa" = c(-2, 1, 1) / 3? Además, cuando los análisis
de ANOVA se compliquen las cosas se pondrán realmente feas.
Sirva el ejemplo de contrasts para resaltar que los
coeficientes de regresión nos permiten hacer contrastes
planeados. Sin embargo, a partir de ahora calcularemos dichos
contrastes con emmeans:
emmeans(v_model, "Species") %>%
# ¡Fíjate que los denominadores ahora sí son comprensibles!
contrast(method = list('V-Setosa' = c(-1, 0.5, 0.5), 'I - II' = c(0, 1, -1)),
infer = c(TRUE, TRUE))
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## V-Setosa 1.256 0.0892 147 1.080 1.432 14.086 <.0001
## I - II -0.652 0.1030 147 -0.855 -0.449 -6.333 <.0001
##
## Confidence level used: 0.95
En general, ANOVA asume:
plot(v_model, which = c(1, 2), ask=FALSE)
# o bien
plot(check_normality(v_model), type = "qq", detrend = TRUE)
check_homogeneity(v_model) # oooooohhhhhhh nooooooooooooo :(
## Warning: Variances differ between groups (Bartlett Test, p = 0.000).
El dataset anxiety proporciona la puntuación de
ansiedad, medida antes y después de aplicar un tratamiento contra la
ansiedad, de tres grupos de personas que practican ejercicios físicos en
diferentes niveles (grp1: basal, grp2: moderado y grp3: alto). Aunque no
tenemos ninguna hipótesis específica, hagamos un análisis de los
datos…
anxiety = read.csv("data/anxiety.csv")
head(anxiety)
## id group pretest posttest
## 1 1 grp1 14.1 14.1
## 2 2 grp1 14.5 14.3
## 3 3 grp1 15.7 14.9
## 4 4 grp1 16.0 15.3
## 5 5 grp1 16.5 15.7
## 6 6 grp1 16.9 16.2
# 1) Vis...
ggplot(anxiety, aes(x = pretest, y = posttest, col = group)) +
geom_point() +
geom_smooth(method="lm")
## `geom_smooth()` using formula 'y ~ x'
# 2) Contrates
# No tenemos hipótesis :(
# 3) Anova + asunciones
anxiety_lm = lm(posttest ~ pretest + group, anxiety)
plot(anxiety_lm, which = c(1, 2), ask=FALSE)
anxiety_aov = Anova(anxiety_lm, type=3)
print(anxiety_aov)
## Anova Table (Type III tests)
##
## Response: posttest
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.285 1 1.6808 0.2021
## pretest 101.289 1 598.3210 <2e-16 ***
## group 74.023 2 218.6293 <2e-16 ***
## Residuals 6.941 41
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
eta_squared(anxiety_aov)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## pretest | 0.94 | [0.90, 1.00]
## group | 0.91 | [0.87, 1.00]
##
## - One-sided CIs: upper bound fixed at (1).
# 4) Posthoc tests!
# ...
¡Ojo, queremos comparar diferencias entre las medias ajustadas! La
función pairwise.t.test() no usará medias ajustadas, por lo que debemos
emplear emmeans:
pairs(
emmeans(anxiety_lm, "group", adjust = "Tukey"),
infer = c(TRUE, TRUE)
)
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## grp1 - grp2 0.641 0.151 41 0.273 1.01 4.235 0.0004
## grp1 - grp3 2.985 0.150 41 2.619 3.35 19.862 <.0001
## grp2 - grp3 2.344 0.151 41 1.976 2.71 15.517 <.0001
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 3 estimates
## P value adjustment: tukey method for comparing a family of 3 estimates
Los diseño de ANOVA factoriales (factorial = más de un factor) permiten el efecto individual y conjunto de uno o más factores. Podemos distinguir varios tipos de análisis factoriales…
… que, como veremos, plantean ciertas diferencias en los modelos. En cualquier caso, en los diseños factoriales lo que primero debemos comprender es el concepto de interacción.
Estudiamos el efecto de tres drogas sobre el tiempo de reacción (una de ellas placebo) teniendo en cuenta además el sexo de los pacientes que toman el medicamento. Supongamos que el efecto de las drogas y edad se mide en términos de reducción del tiempo de reacción a algún estímulo y que se obtienen los resultados del fichero “drugs_1.csv”. Visualiza el efecto de las drogas y sexo en los tiempos de reacción y propón un modelo.
drugs_df_1 =
read.csv("data/drugs_1.csv") %>%
mutate(drug = factor(drug), sex = factor(sex))
# interaction.plot(drugs_df$sex, drugs_df$drug, response = drugs_df$response_time)
ggplot(drugs_df_1, aes(x=sex, y=response_time, col=drug)) +
stat_summary() +
stat_summary(geom='line', aes(group=drug))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
drugs_model_1 = lm(response_time ~ sex + drug, data = drugs_df_1)
drugs_df_1$predictions = predict(drugs_model_1)
ggplot(drugs_df_1, aes(x=sex, y=predictions, col=drug)) +
stat_summary() +
stat_summary(geom='line', aes(group=drug))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
Repite el ejercicio anterior para los datos experimentales “drugs_2.csv”.
drugs_df_2 =
read.csv("data/drugs_2.csv") %>%
mutate(drug = factor(drug), sex = factor(sex))
# interaction.plot(drugs_df$sex, drugs_df$drug, response = drugs_df$response_time)
ggplot(drugs_df_2, aes(x=sex, y=response_time, col=drug)) +
stat_summary() +
stat_summary(geom='line', aes(group=drug))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
drugs_model_2 = lm(response_time ~ sex + drug, data = drugs_df_2)
drugs_df_2$predictions = predict(drugs_model_2)
ggplot(drugs_df_2, aes(x=sex, y=predictions, col=drug)) +
stat_summary() +
stat_summary(geom='line', aes(group=drug))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
Uuuuups, las predicciones son malíiiiiiisimas… El modelo no es capaz de capturar las interacciones presentes en los datos.
Para modelar interacciones…
drugs_model_2 = lm(response_time ~ sex * drug, data = drugs_df_2)
# o de forma equivalente
#drugs_model_2 = lm(response_time ~ sex + drug + sex:drug, data = drugs_df_2)
drugs_df_2$predictions = predict(drugs_model_2)
ggplot(drugs_df_2, aes(x=sex, y=predictions, col=drug)) +
stat_summary() +
stat_summary(geom='line', aes(group=drug))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
Una vez aclarado el concepto de interacción podemos completar el
análisis para drugs_df_2.
# 2) Contrastes: aunque hemos dejado contrasts de lado, es recomendable
# planear los contrastes antes de correr ANOVA.
# ¡Ojo! ahora creamos listas de contrastes
drugs_contrasts = list(
"_Drugs-Placebo" = c(1, 1, -2),
"_A - B" = c(1, -1, 0)
)
# 3) ANOVA
drugs_model_2 = lm(response_time ~ sex * drug, data = drugs_df_2)
print(Anova(drugs_model_2, type = 3))
## Anova Table (Type III tests)
##
## Response: response_time
## Sum Sq Df F value Pr(>F)
## (Intercept) 2800.49 1 3985.48 < 2.2e-16 ***
## sex 609.76 1 867.77 < 2.2e-16 ***
## drug 1091.94 2 776.99 < 2.2e-16 ***
## sex:drug 660.92 2 470.29 < 2.2e-16 ***
## Residuals 16.86 24
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
eta_squared(drugs_model_2)
## # Effect Size for ANOVA (Type I)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------
## sex | 0.71 | [0.53, 1.00]
## drug | 0.99 | [0.99, 1.00]
## sex:drug | 0.98 | [0.96, 1.00]
##
## - One-sided CIs: upper bound fixed at (1).
# Omitimos los plot de comprobación por sencillez...(comprueba que son correctos)
# plot(drugs_model_2, ask=FALSE)
# 4) contrastes
# Contrastes principales para drogas
emmeans(drugs_model_2, ~ drug) %>%
contrast(method = list("drugs" = drugs_contrasts))
## NOTE: Results may be misleading due to involvement in interactions
## contrast estimate SE df t.ratio p.value
## drugs._Drugs-Placebo 40.132 0.649 24 61.808 <.0001
## drugs._A - B -0.997 0.375 24 -2.659 0.0137
##
## Results are averaged over the levels of: sex
# Contrastes para interacciones
emmeans(drugs_model_2, ~ drug | sex) %>%
contrast(interaction = list("drugs" = drugs_contrasts, "sex" = "consec"), by=NULL)
## drug_custom sex_consec estimate SE df t.ratio p.value
## _Drugs-Placebo Male - Female 18.9 1.30 24 14.588 <.0001
## _A - B Male - Female 20.2 0.75 24 26.977 <.0001
Lo más interesante del código anterior es reflexionar sobre cada uno
de los contrastes y su significado. Los contrastes básicos
(sex, Drugs-Placebo, A - B)
deberían ser claros pero, ¿qué significan los contrastes para
interacciones?
sex1:drug_Drugs-Placebo: El efecto
Drugs-Placebo es diferente en hombres y mujeres? Es decir,
¿el efecto de placebo Vs drogas es comparable en hombres y mujeres?sex1:drug_A - B: El efecto drug_A - B es
diferente en hombres y mujeres? Es decir, ¿el efecto droga A Vs droga B
es comparable en hombres y mujeres?Veamos cómo visualizar que las interacciones entre contrastes y sexo son significativas…
contraste_1 = drugs_df_2 %>% mutate(drug = fct_recode(drug, 'Drug' = 'A', 'Drug' = 'B'))
contraste_2 = drugs_df_2 %>% filter(drug != "Placebo") %>% mutate(drug = fct_drop(drug))
ggplot(contraste_1, aes(x=drug, y=response_time, col=sex)) +
stat_summary() +
stat_summary(geom='line', aes(group=sex))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
ggplot(contraste_2, aes(x=drug, y=response_time, col=sex)) +
stat_summary() +
stat_summary(geom='line', aes(group=sex))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
En ambos casos, las pendientes de los grupos son distintas, lo que apoya que hay interacciones en los datos.
Fíjate que, si las interacciones son significativas, la interpretación de los efectos principales no tiene sentido. Por ejemplo, en el contraste 2, la droga B aumenta el tiempo de respuesta para las mujeres, pero lo disminuye para hombres. Por tanto, responder a ¿disminuye la droga B el tiempo de respuesta (para hombres y mujeres)? da una imagen incompleta del problema.
NO INTERPRETES LOS EFECTOS PRINCIPALES SI LA INTERACCIÓN ES SIGNIFICATIVA.
ggplot(drugs_df_2, aes(x=sex, y=response_time, col = drug)) +
stat_summary() +
stat_summary(geom="line", aes(group = drug))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
drugs_emms = emmeans(drugs_model_2, ~ drug | sex)
contrast(drugs_emms, interaction = list(drug="pairwise", sex='consec'), by=NULL,
adjust = "fdr")
## drug_pairwise sex_consec estimate SE df t.ratio p.value
## A - B Male - Female 20.226 0.75 24 26.977 <.0001
## A - Placebo Male - Female 19.585 0.75 24 26.122 <.0001
## B - Placebo Male - Female -0.641 0.75 24 -0.855 0.4012
##
## P value adjustment: fdr method for 3 tests
Considera un experimento en el que a cada paciente se le realizan tres mediciones: una al comenzar la prueba, otra medición tras probar el tratamiento A, y otra medición tras probar el tratamiento B…
Los tests F en ANOVA asumen que las medidas son independientes. Cuando se utilizan medidas repetidas como en el ejemplo anterior, se viola esta suposición. Por tanto, cuando tenemos ANOVA con medidas repetidas debemos “cambiar” este supuesto por otro válido. Se asume que la relación entre pares de condiciones experimentales es similar. A este supuesto se le conoce como esfericidad. De forma práctica…
Desafortunadamente, el test car::Anova no incorpora
correcciones contra la violación de esfericidad por lo que lo
cambiaremos por afex::aov_4.
Existe evidencia de que las actitudes hacia ciertos estímulos se pueden cambiar utilizando imágenes positivas. Como parte de una iniciativa para detener el consumo excesivo de alcohol en los adolescentes, el gobierno financió un estudio para determinar si imágenes negativas podrían usarse para hacer que las actitudes de los adolescentes hacia el alcohol sean más negativas. Cada participante en cada estudio vio un total de nueve anuncios simulados en tres sesiones. En una sesión, cada partipante ve tres anuncios (de cerveza, vino o agua); cada anuncio tiene una actitud asociada (positiva, negativa o neutra). Al cabo de las 3 sesiones, las variables se cruzan completamente produciendo 9 condiciones experimentales. Después de cada anuncio, se pidió a los participantes que calificaran las bebidas en una escala que iba desde -100 (me disgusta mucho) pasando por 0 (neutral) hasta 100 (me gusta mucho). El orden de los anuncios fue aleatorio, al igual que el orden en que las personas participaron en las tres sesiones. Datos en “data/attitude_long.dat” (o “attitude.dat” para un reto).
# Los datos en attitude.dat no son adecuados para anova... hay que arreglarlo!
attitude_df = read.table("data/attitude.dat", header = TRUE)
attitude_df_long =
attitude_df %>%
pivot_longer(beer_pos:water_neu, names_to = 'name', values_to = 'attitude') %>%
separate(name, into=c('drink', 'imagery'), sep = '_') %>%
mutate_if(is.character, as.factor)
# 1) plot
attitude_df_long = read.table("data/attitude_long.dat", header=TRUE,
stringsAsFactors = TRUE)
head(attitude_df_long)
## participant drink imagery attitude
## 1 P1 beer pos 1
## 2 P1 beer neg 6
## 3 P1 beer neu 5
## 4 P1 wine pos 38
## 5 P1 wine neg -5
## 6 P1 wine neu 4
ggplot(attitude_df_long, aes(x=drink, y=attitude, fill=imagery)) + geom_boxplot()
ggplot(attitude_df_long, aes(x=drink, y=attitude, col=imagery)) +
stat_summary() +
stat_summary(geom="line", aes(group=imagery))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
library('afex')
## Loading required package: lme4
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
## ************
## Welcome to afex. For support visit: http://afex.singmann.science/
## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - NEWS: emmeans() for ANOVA models now uses model = 'multivariate' as default.
## - Get and set global package options with: afex_options()
## - Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************
##
## Attaching package: 'afex'
## The following object is masked from 'package:lme4':
##
## lmer
# 2) Contrastes: usaremos agua para comparar las bebidas alcohólicas y no-alcohólicas
drinks_contrasts = list(
"Alcohol-Water" = c(1, -2, 1),
"Beer-Wine" = c(1, 0, -1)
)
imagery_contrasts = list(
"Others-Neutral" = c(1, -2, 1),
"Pos-Neg" = c(-1, 0, 11)
)
# 3) ANOVA y asunciones
attitude_model = aov_4(attitude ~ drink * imagery + (drink * imagery | participant),
data = attitude_df_long, observed=NULL)
# El effect size se corresponde con ges (generalized eta squared)
summary(attitude_model)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 11218.0 1 1920.1 19 111.005 2.255e-09 ***
## drink 2092.3 2 7785.9 38 5.106 0.01086 *
## imagery 21628.7 2 3352.9 38 122.565 < 2.2e-16 ***
## drink:imagery 2624.4 4 2906.7 76 17.155 4.589e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## drink 0.26724 0.00001
## imagery 0.66210 0.02445
## drink:imagery 0.59504 0.43566
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## drink 0.57711 0.02977 *
## imagery 0.74744 1.757e-13 ***
## drink:imagery 0.79840 1.900e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## drink 0.5907442 2.881391e-02
## imagery 0.7968420 3.142804e-14
## drink:imagery 0.9785878 6.809640e-10
# asunciones
plot(check_normality(attitude_model), "qq", detrend = TRUE)
check_sphericity(attitude_model) # ya lo sabíamos :(
## Warning: Sphericity violated for:
## - drink (p < .001)
## - imagery (p = 0.024).
# 4.a) Contrastes ...
emm = emmeans(attitude_model, ~ drink | imagery)
interaction_contrasts = list(
"drink" = drinks_contrasts,
"imagery" = imagery_contrasts
)
# 4.b) Test post-hoc
contrast(emm, interaction = interaction_contrasts, by = NULL)
## drink_custom imagery_custom estimate SE df t.ratio p.value
## Alcohol-Water Others-Neutral -11.4 7.53 19 -1.521 0.1446
## Beer-Wine Others-Neutral 15.4 4.77 19 3.242 0.0043
## Alcohol-Water Pos-Neg 116.8 53.33 19 2.189 0.0413
## Beer-Wine Pos-Neg -63.8 36.02 19 -1.770 0.0928
attitude_emms = emmeans(attitude_model, ~ imagery | drink)
contrast(attitude_emms, interaction = list(imagery = "pairwise", drink = "consec"),
by=NULL, adjust = "fdr")
## imagery_pairwise drink_consec estimate SE df t.ratio p.value
## neg - neu water - beer -6.00 2.31 19 -2.599 0.0265
## neg - pos water - beer -10.00 3.26 19 -3.063 0.0128
## neu - pos water - beer -4.00 3.19 19 -1.255 0.2695
## neg - neu wine - water -12.10 2.33 19 -5.187 0.0003
## neg - pos wine - water -10.75 2.65 19 -4.056 0.0020
## neu - pos wine - water 1.35 2.78 19 0.485 0.6334
##
## P value adjustment: fdr method for 6 tests
# Visualización de apoyo a los contrastes/tests
ggplot(attitude_df_long, aes(x = drink, y = attitude, col = imagery)) +
stat_summary() +
stat_summary(geom="line", aes(group = imagery))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
Como nota para cerrar ANOVA, una alternativa más flexible a estos
modelos son los Linear Mixture Models (LMMs), que pueden manejar datos
no-independientes, no-esféricos, etc. Los paquetes principales en
R son nlme y lme4.
Como parte de una investigación de agentes tóxicos, se asignaron 48 ratas a 3 venenos (I,II,III) y 4 tratamientos (A,B,C,D). La respuesta fue el tiempo de supervivencia en decenas de horas. Realiza una análisis ANOVA completo de los datos.
Hemos visto que la Normalidad de los datos es una asunción importante de los métodos vistos hasta el momento. ¿Qué hacer si los datos violan claramente la normalidad?
Los metodos no paramétricos permiten realizar análisis estadísticos sin asumir ninguna distribución de los datos. A cambio, pierden potencia en los contrastes.
A grandes rasgos tenemos… * Como alternativa a t-test independiente
\(\rightarrow\) Mann–Whitney test o
Wilcoxon’s rank-sum test \(\rightarrow\) wilcox.test. *
Como alternativa a t-test apareado \(\rightarrow\) Wilcoxon signed-rank test
\(\rightarrow\)
wilcox.test. * Como alternativa a One-way ANOVA \(\rightarrow\) Kruskal–Wallis test \(\rightarrow\) kruskal.test. *
Como alternativa a Repeated Measures ANOVA \(\rightarrow\) Friedman’s ANOVA \(\rightarrow\)
friedman.test.
El gasto cardı́aco (litros/minuto) se midió por termodilución en una muestra aleatoria simple de 15 pacientes operados de corazón. Deseamos saber si podemos concluir que la media de la población es diferente de 5.05.
data = c(4.91, 4.10, 6.74, 7.27, 7.42, 7.50, 6.56, 4.64, 5.98, 3.14, 3.23,
5.80,6.17,5.39, 5.77)
w_test = wilcox.test(data, mu = 5.05)
print(w_test)
##
## Wilcoxon signed rank exact test
##
## data: data
## V = 86, p-value = 0.1514
## alternative hypothesis: true location is not equal to 5.05
w_test %>% rank_biserial()
## r (rank biserial) | 95% CI
## ---------------------------------
## 0.43 | [-0.11, 0.78]
##
## - Deviation from a difference of 5.05.
En un estudio, se midieron los estreses hemodinámicos en sujetos sometidos a colecistectomı́a laparoscópica. Una variable de interés es el volumen diastólico final ventricular (LVEDV) medido en mililitros. Los datos se encuentra en “lvedv.tsv”. El “baseline” se refiere a una medición tomada 5 minutos después de la inducción de la anestesia, y el término “5 minutes” se refiere a una medición tomada 5 minutos después del “baseline”. ¿Hay cambios en los niveles de LVEDV durante la intervención?
lvedv = read.table("data/lvedv.tsv", header = TRUE)
w_test = wilcox.test(lvedv$baseline, lvedv$X5_minutes, paired = TRUE)
## Warning in wilcox.test.default(lvedv$baseline, lvedv$X5_minutes, paired = TRUE):
## cannot compute exact p-value with ties
print(w_test)
##
## Wilcoxon signed rank test with continuity correction
##
## data: lvedv$baseline and lvedv$X5_minutes
## V = 11.5, p-value = 0.1139
## alternative hypothesis: true location shift is not equal to 0
w_test %>% rank_biserial()
## r (rank biserial) | 95% CI
## ---------------------------------
## -0.58 | [-0.88, 0.03]
Se ha diseñado un experimento para evaluar los efectos de la inhalación prolongada de óxido de cadmio. Quince animales de laboratorio sirvieron como sujetos experimentales, mientras que 10 animales similares sirvieron como controles. La variable de interés fue el nivel de hemoglobina después del experimento. Los resultados se encuentran en “cadmium.csv”. Deseamos saber si podemos concluir que la inhalación prolongada de óxido de cadmio reduce el nivel de hemoglobina
dat = read.table("data/cadmium.csv", header = TRUE, sep = ",")
w_test = wilcox.test(
dat %>% filter(class == "exposed") %>% .$hemoglobin,
dat %>% filter(class == "unexposed") %>% .$hemoglobin,
alternative = "less"
)
## Warning in wilcox.test.default(dat %>% filter(class == "exposed") %>% .
## $hemoglobin, : cannot compute exact p-value with ties
print(w_test)
##
## Wilcoxon rank sum test with continuity correction
##
## data: dat %>% filter(class == "exposed") %>% .$hemoglobin and dat %>% filter(class == "unexposed") %>% .$hemoglobin
## W = 25, p-value = 0.003004
## alternative hypothesis: true location shift is less than 0
w_test %>% rank_biserial()
## r (rank biserial) | 95% CI
## ----------------------------------
## -0.67 | [-1.00, -0.39]
##
## - One-sided CIs: lower bound fixed at (-1).
En PlantGrowth tenemos los resultados de un experimento
para comparar los rendimientos (medidos por el peso en seco de las
plantas) obtenidos bajo un control y dos tratamientos diferentes para
hacer crecer a las plantas.
data("PlantGrowth")
head(PlantGrowth)
## weight group
## 1 4.17 ctrl
## 2 5.58 ctrl
## 3 5.18 ctrl
## 4 6.11 ctrl
## 5 4.50 ctrl
## 6 4.61 ctrl
# Alternativa: kruskal.test(weight ~ group, data = PlantGrowth)
kw_test = kruskal.test(PlantGrowth$weight, PlantGrowth$group)
kw_test %>% rank_epsilon_squared()
## Epsilon2 (rank) | 95% CI
## ------------------------------
## 0.28 | [0.11, 1.00]
##
## - One-sided CIs: upper bound fixed at (1).
pairwise.wilcox.test(PlantGrowth$weight, PlantGrowth$group,
p.adjust.method = "BH")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: PlantGrowth$weight and PlantGrowth$group
##
## ctrl trt1
## trt1 0.199 -
## trt2 0.095 0.027
##
## P value adjustment method: BH
El dataset selfesteem contiene la puntuación de
autoestima de 10 personas con obesidad en tres momentos durante una
dieta. El objetivo es determinar si su autoestima mejoró con la
dieta.
selfesteem = read.csv("data/selfesteem.csv")
head(selfesteem)
## id time score
## 1 1 t1 4.005027
## 2 1 t2 5.182286
## 3 1 t3 7.107831
## 4 2 t1 2.558124
## 5 2 t2 6.912915
## 6 2 t3 6.308434
ggplot(selfesteem, aes(x=time, y=score)) + geom_boxplot()
# Alternativa: friedman.test(score ~ time | id, selfesteem)
f_test = friedman.test(selfesteem$score, groups = selfesteem$time, blocks = selfesteem$id)
f_test %>% kendalls_w()
## Kendall's W | 95% CI
## --------------------------
## 0.91 | [0.79, 1.00]
##
## - One-sided CIs: upper bound fixed at (1).
pairwise.wilcox.test(selfesteem$score, selfesteem$time, paired = TRUE,
p.adjust.method = "bonferroni")
##
## Pairwise comparisons using Wilcoxon signed rank exact test
##
## data: selfesteem$score and selfesteem$time
##
## t1 t2
## t2 0.0059 -
## t3 0.0059 0.0117
##
## P value adjustment method: bonferroni
La carpeta “meta” contiene los datos de varios papers de los que se
desea hacer un meta-análisis. Como primer paso, deseas verificar que los
resultados que obtienes sobre los datos son los mismos que se han
reportado en los papers originales. Para ello, leerás los datos (todos
tienen dos columnas: response y group), ejecutarás un T-test y generarás
un gráfico donde se comparen las respuestas para los niveles de group.
Los gráficos resultantes se guardarán en la carpeta “meta/images” y los
p-valores de los tests en un solo data.frame
p-values.
base_path = "data/meta"
images_path = file.path(base_path, "images")
print(images_path)
## [1] "data/meta/images"
dir.create(images_path, showWarnings = FALSE)
files = list.files(base_path, pattern=".csv")
p_values = data.frame()
for (file in files) {
df = read.csv(file.path(base_path, file))
test = t.test(response ~ group, df)
result = data.frame("filename" = file, "p_value" = test$p.value)
p_values = rbind(p_values, result)
plot = ggplot(df, aes(x = group, y = response, fill = group)) + geom_boxplot()
image_name = gsub(".csv", ".png", file)
ggsave(file.path(images_path, image_name), plot)
}
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
print(p_values)
## filename p_value
## 1 paper-1.csv 2.050053e-01
## 2 paper-10.csv 1.170465e-01
## 3 paper-11.csv 3.461138e-07
## 4 paper-12.csv 8.363042e-02
## 5 paper-13.csv 3.838954e-01
## 6 paper-14.csv 2.498351e-05
## 7 paper-15.csv 5.605239e-02
## 8 paper-16.csv 1.036900e-05
## 9 paper-17.csv 6.960640e-01
## 10 paper-18.csv 1.825461e-13
## 11 paper-19.csv 1.300555e-07
## 12 paper-2.csv 5.266617e-01
## 13 paper-20.csv 6.377956e-02
## 14 paper-3.csv 1.030685e-05
## 15 paper-4.csv 2.305727e-04
## 16 paper-5.csv 1.445790e-07
## 17 paper-6.csv 1.603187e-01
## 18 paper-7.csv 1.116642e-09
## 19 paper-8.csv 1.279665e-14
## 20 paper-9.csv 1.254387e-07
write.csv(p_values, "p_values.csv", row.names = FALSE)
J. Després, D. Homme, M. Pouliot, A. Tremblay, and C. Bouchard, ``Estimation of Deep Abdominal Adipose-Tissue Accumulation from Simple Anthropometric Measurements in Men’’ American Journal of Clinical Nutrition, 54 (1991), 471–477.↩︎